00010	BEGIN "FAST"
00020	COMMENT		This program does a series of fast Fourier transforms ON
00030			Segments of  data taken from an
00040			input file specified from  the terminal 
00050			and creates an output file containing  a series of
00060			power spectra. The size of the segments to be used
00070			used is specified by typing in a value for M,
00080			which must be ≤ 8;
00090	
00100	REQUIRE "DPYSUB.HDR[1,PDQ]" SOURCE_FILE;
00110	DEFINE DPYSIZ="1000";
00120	INTEGER ARRAY DPYBUF[1:DPYSIZ];
00130	
00140	EXTERNAL STRING PROCEDURE DATIME;
00150	EXTERNAL PROCEDURE TIMSET;
00160	EXTERNAL STRING PROCEDURE STRIN(STRING S);
00170	EXTERNAL INTEGER PROCEDURE ININT(STRING S);
00180	EXTERNAL REAL PROCEDURE RUNTIM;
00190	EXTERNAL STRING PROCEDURE INCHWL;
00200	EXTERNAL PROCEDURE WAIT(INTEGER I);
00210	
00220	FORTRAN REAL PROCEDURE ALOG10(REAL X);
00230	FORTRAN REAL PROCEDURE COS(REAL X);
00240	FORTRAN REAL PROCEDURE SIN(REAL X);
00250	FORTRAN REAL PROCEDURE SQRT(REAL X);
00260	
00270	REAL ARRAY A,B,C[0:256];
00280	REAL ARRAY W[0:128];
00290	
00300	STRING FILEI,TFILEI,CHAR1,CHAR2,CHAR3,CHAR4,CHAR5,CHAR6,CHAR7,LOGFRE,CEPS;
00310	
00320	INTEGER ARRAY LFILE[0:'177];
00330	INTEGER ARRAY D[0:256];
00340	INTEGER ARRAY DATBUF[0:1023];
00350	
00360	INTEGER AIFORM,AOFORM,EOF,IEOF,I,J,K,M,TM,N,BPT,SEGC,SEGCNT,SEGTOT,SAMINT;
00370	INTEGER XPOS,YPOS,MAX;
00380	EXTERNAL REAL RTIM,ATIM;
00390	
00400	DEFINE DATSIZ="1024";
00410	DEFINE BPS="12";
00420	DEFINE LEFT="24";
00430	DEFINE DI="2↑24";
00440	
00450	LABEL START;
00460	
     

00010	PROCEDURE UNPACK(REAL ARRAY A,B;INTEGER ARRAY DATBUF;REFERENCE INTEGER BPT;INTEGER N,SAMINT);
00020	BEGIN
00030	COMMENT		This procedure accepts a segment of digitized speech of
00040			length 2↑(m+1). It expands this segment and
00050			loads it into A andB in real form;
00060		INTEGER I,J;
00070		FOR I←0 STEP 1 UNTIL N-1 DO
00080		  BEGIN
00090		    A[I]←((ILDB(BPT) LSH LEFT)%DI);
00100		    IF (SAMINT)=2 THEN J←ILDB(BPT);
00110		  END;
00120		FOR I←0 STEP 1 UNTIL N-1 DO
00130		  BEGIN
00140		    B[I]←((ILDB(BPT) LSH LEFT)%DI);
00150		    IF (SAMINT)=2 THEN J←ILDB(BPT);
00160		  END;
00170	END "UNPACK";
00180	
00190	PROCEDURE DPYDMP;
00200	BEGIN	STRING FILE;INTEGER DSIZ,CCCHN;
00210		IF ¬(FILE←STRIN("CAL-COMP FILE=")) THEN RETURN;
00220		OPEN(CCCHN←GETCHAN,"DSK",'14,0,1,0,0,0);
00230		ENTER(CCCHN,FILE,0);
00240		DPYPARS;DSIZ←DPYBUF[2]+3;
00250		ARRYOUT(CCCHN,DPYBUF[1],DSIZ);
00260		RELEASE(CCCHN);
00270	END "DPYDMP";
00280	
00290	
00300	PROCEDURE FFT(REAL ARRAY A,B;INTEGER N,M,KS);
00310	BEGIN
00320	COMMENT COMPUTES THE FFT FOR ONE VARIABLE OF DIMENSION 2↑M;
00330	INTEGER K0,K1,K2,K3,SPAN,J,JJ,K,KB,KN,MM,MK;
00340	REAL RAD,C1,C2,C3,S1,S2,S3,CK,SK,SQ;
00350	REAL A0,A1,A2,A3,B0,B1,B2,B3;
00360	INTEGER ARRAY C[0:M];
00370	LABEL L,L2,L3,L4,L5,L6;
00380	SQ←0.707106781187;
00390	SK←0.382683432366;
00400	CK←0.92387953251;
00410	C[M]←KS; MM←(M%2)*2; KN←0;
00420	FOR K←M-1 STEP -1 UNTIL 0 DO C[K]←C[K+1]/2;
00430	RAD←6.28318530718/(C[0]*KS); MK←M-5;
00440	L:  KB←KN; KN←KN+KS;
00450	IF MM≠N THEN
00460	BEGIN
00470		K2←KN;  K0←C[MM]+KB;
00480	L2:	K2←K2-1; K0←K0-1;
00490		A0←A[K2]; B0←B[K2];
00500		A[K2]←A[K0]-A0; A[K0]←A[K0]+A0;
00510		B[K2]←B[K0]-B0; B[K0]←B[K0]+B0;
00520		IF K0>KB THEN GO TO L2;
00530	END;
00540	C1←1.0; S1←0;
00550	JJ←0; K←MM-2; J←3;
00560	IF K≥0 THEN GO TO L4 ELSE GO TO L6;
00570	L3:  IF C[J]≤JJ THEN
00580	BEGIN
00590		JJ←JJ-C[J]; J←J-1;
00600		IF C[J]≤JJ THEN
00610		BEGIN
00620			JJ←JJ-C[J]; J←J-1; K←K+2;
00630			GO TO L3;
00640		END
00650	END;
00660	JJ←C[J]+JJ; J←3;
00670	L4:  SPAN←C[K];
00680	IF JJ≠0 THEN
00690	BEGIN
00700		C2←JJ*SPAN*RAD; C1←COS(C2); S1←SIN(C2);
00710	L5:	C2←C1↑2-S1↑2; S2←2.0*C1*S1;
00720		C3←C2*C1-S2*S1; S3←C2*S1+S2*C1;
00730	END;
00740	FOR K0←KB+SPAN-1 STEP -1 UNTIL KB DO
00750	BEGIN
00760		K1←K0+SPAN; K2←K1+SPAN; K3←K2+SPAN;
00770		A0←A[K0]; B0←B[K0];
00780		IF S1=0 THEN
00790		BEGIN
00800		A1←A[K1]; B1←B[K1];
00810		A2←A[K2]; B2←B[K2];
00820		A3←A[K3]; B3←B[K3];
00830		END
00840		ELSE
00850		BEGIN
00860		A1←A[K1]*C1-B[K1]*S1;
00870		B1←A[K1]*S1+B[K1]*C1;
00880		A2←A[K2]*C2-B[K2]*S2;
00890		B2←A[K2]*S2+B[K2]*C2;
00900		A3←A[K3]*C3-B[K3]*S3;
00910		B3←A[K3]*S3+B[K3]*C3;
00920		END;
00930		A[K0]←A0+A2+A1+A3; B[K0]←B0+B2+B1+B3;
00940		A[K1]←A0+A2-A1-A3; B[K1]←B0+B2-B1-B3;
00950		A[K2]←A0-A2-B1+B3; B[K2]←B0-B2+A1-A3;
00960		A[K3]←A0-A2+B1-B3; B[K3]←B0-B2-A1+A3;
00970	END;
00980	IF K>0 THEN BEGIN K←K-2;  GO TO L4; END;
00990	KB←K3+SPAN;
01000	IF KB<KN THEN
01010	BEGIN
01020		IF J=0 THEN BEGIN K←2; J←MK; GO TO L3; END;
01030		J←J-1; C2←C1;
01040		IF J=1 THEN
01050		BEGIN C1←C1*CK+S1*SK; S1←S1*CK-C2*SK; END
01060		ELSE BEGIN C1←(C1-S1)*SQ; S1←(C2+S1)*SQ; END;
01070		GO TO L5;
01080	END;
01090	L6: IF KN<N THEN GO TO L;
01100	END "FFT";
01110	
01120	
01130	
01140	
01150	PROCEDURE REVFFT(REAL ARRAY A,B;INTEGER N,M,KS);
01160	BEGIN
01170	COMMENT COMPUTES THE FFT FOR ONE VARIABLE OF DIMENSION 2↑M IN A
01180		MULTIVARIATE TRANSFORM.
01190		IF N=2↑M AND K=1 THEN A SINGLE-VARIATE TRANSFORM IS COMPUTED;
01200	INTEGER K0,K1,K2,K3,K4,SPAN,NN,J,JJ,K,KB,NT,KN,MK;
01210	REAL RAD,C1,C2,C3,S1,S2,S3,CK,SK,SQ;
01220	REAL A0,A1,A2,A3,B0,B1,B2,B3,RE,IM;
01230	INTEGER ARRAY C[0:M];
01240	LABEL L,L2,L3,L4,L5,L6;
01250	SQ←0.707106781187;
01260	SK←0.382683432366;
01270	CK←0.92387953251;
01280	C[0]←KS; KN←0; K4←4*KS; MK←M-4;
01290	FOR K←1 STEP 1 UNTIL M DO C[K]←KS←KS+KS;
01300	RAD←3.1415926536/(C[0]*KS);
01310	L:  KB←KN+K4; KN←KN+KS;
01320	IF M=1 THEN GO TO L5;
01330	K←JJ←0; J←MK; NT←3;
01340	C1←1.0; S1←0;
01350	L2:  SPAN←C[K];
01360	IF JJ≠0 THEN
01370	BEGIN
01380		C2←JJ*SPAN*RAD; C1←COS(C2); S1←SIN(C2);
01390	L3:	C2←C1↑2-S1↑2; S2←2*C1*S1;
01400		C3←C2*C1-S2*S1; S3←C2*S1+S2*C1;
01410	END 
01420	ELSE S1←0;
01430	K3←KB-SPAN;
01440	L4:  K2←K3-SPAN; K1←K2-SPAN; K0←K1-SPAN;
01450		A0←A[K0]; B0←B[K0];
01460		A1←A[K1]; B1←B[K1];
01470		A2←A[K2]; B2←B[K2];
01480		A3←A[K3]; B3←B[K3];
01490		A[K0]←A0+A1+A2+A3; B[K0]←B0+B1+B2+B3;
01500	IF S1=0 THEN
01510	BEGIN
01520		A[K1]←A0-A1-B2+B3; B[K1]←B0-B1+A2-A3;
01530		A[K2]←A0+A1-A2-A3; B[K2]←B0+B1-B2-B3;
01540		A[K3]←A0-A1+B2-B3; B[K3]←B0-B1-A2+A3;
01550	END
01560	ELSE
01570	BEGIN
01580		RE←A0-A1-B2+B3; IM←B0-B1+A2-A3;
01590		A[K1]←RE*C1-IM*S1; B[K1]←RE*S1+IM*C1;
01600		RE←A0+A1-A2-A3; IM←B0+B1-B2-B3;
01610		A[K2]←RE*C2-IM*S2;  B[K2]←RE*S2+IM*C2;
01620		RE←A0-A1+B2-B3; IM←B0-B1-A2+A3;
01630		A[K3]←RE*C3-IM*S3;  B[K3]←RE*S3+IM*C3;
01640	END;
01650	K3←K3+1; IF K3<KB THEN GO TO L4;
01660	NT←NT-1;
01670	IF NT≥0 THEN
01680	BEGIN
01690		C2←C1;
01700		IF NT=1 THEN
01710		BEGIN C1←C1*CK+S1*SK; S1←S1*CK-C2*SK; END
01720		ELSE BEGIN C1←(C1-S1)*SQ; S1←(C2+S1)*SQ; END;
01730		KB←KB+K4; IF KB≤KN THEN GO TO L3 ELSE GO TO L5;
01740	END;
01750	IF NT=-1 THEN BEGIN K←2; GO TO L2; END;
01760	IF C[J]≤JJ THEN
01770	BEGIN
01780		JJ←JJ-C[J]; J←J-1;
01790		IF C[J]≤JJ THEN
01800		BEGIN JJ←JJ-C[J]; J←J-1; K←K+2; END
01810		ELSE BEGIN JJ←C[J]+JJ; J←MK;END;
01820	END
01830	ELSE BEGIN JJ←C[J]+JJ; J←MK; END;
01840	IF J<MK THEN GO TO L2; K←0; NT←3;
01850	KB←KB+K4;  IF KB≤KN THEN GO TO L2;
01860	L5:  K←(M%2)*2;
01870	IF K≠M THEN
01880	BEGIN
01890		K2←KN; K0←J←KN-C[K];
01900	L6:	K2←K2-1; K0←K0-1;
01910		A0←A[K2]; B0←B[K2];
01920		A[K2]←A[K0]-A0; A[K0]←A[K0]+A0;
01930		B[K2]←B[K0]-B0; B[K0]←B[K0]+B0;
01940		IF K2>J THEN GO TO L6;
01950	END;
01960	IF KN<N THEN GO TO L;
01970	END "REVFFT";
01980	
01990	PROCEDURE REORDER(REAL ARRAY A,B;INTEGER N,M,KS,REEL);
02000	BEGIN
02010	COMMENT PERMUTES DATA FROM NORMAL TO REVERSE BINARY ORDER AND BACK;
02020	INTEGER I,J,JJ,K,KK,KB,K2,KU,LIM,P;
02030	REAL T;
02040	INTEGER ARRAY C,LST[0:M];
02050	LABEL L,L2,L3,L4;
02060	C[M]←KS;
02070	FOR K←M STEP -1 UNTIL 1 DO C[K-1]←C[K]%2;
02080	 P←J←M-1; I←KB←0;
02090	IF REEL THEN
02100	BEGIN
02110		KU←N-2;
02120		FOR K←0 STEP 2 UNTIL KU DO
02130		BEGIN T←A[K+1]; A[K+1]←B[K]; B[K]←T; END;
02140	END
02150	ELSE M←M-1;
02160	LIM←(M+2)%2; IF P≤0 THEN GO TO L4;
02170	L:	KU←K2←C[J]+KB; JJ←C[M-J]; KK←KB+JJ;
02180	L2:	K←KK+JJ;
02190	L3:	T←A[KK]; A[KK]←A[K2]; A[K2]←T;
02200	T←B[KK]; B[KK]←B[K2]; B[K2]←T;
02210	KK←KK+1; K2←K2+1;
02220	IF KK<K THEN GO TO L3;
02230	KK←KK+JJ; K2←K2+JJ;
02240	IF KK<KU THEN GO TO L2;
02250	IF J>LIM THEN
02260	BEGIN
02270		J←J-1; I←I+1;
02280		LST[I]←J; GO TO L;
02290	END;
02300	KB←K2;
02310	IF I>0 THEN
02320	BEGIN J←LST[I]; I←I-1; GO TO L; END;
02330	IF KB<N THEN BEGIN J←P; GO TO L; END;
02340	L4:   ;
02350	END "REORDER";
02360	
02370	
02380	PROCEDURE RTRAN(REAL ARRAY A,B;INTEGER N,EVALUATE);
02390	BEGIN
02400	COMMENT IF EVALUATE IS FALSE THIS PROCEDURE UNSCRAMBLES  THE SINGLE VARIATE
02410	COMPLEX TRANSFORM ;
02420	INTEGER K,NK,NH;
02430	REAL AA,AB,BA,BB,RE,IM,CK,SK,DC,DS,R;
02440	NH←N%2;  R←3.1415926536/N;
02450	DS←SIN(R); R←-(2*SIN(0.5*R))↑2;
02460	DC←-0.5*R; CK←1.0;  SK←0;
02470	IF EVALUATE THEN
02480	BEGIN
02490	CK←-1.0; DC←-DC;
02500	END
02510	ELSE
02520	BEGIN
02530	A[N]←A[0]; B[N]←B[0];
02540	END;
02550	FOR K←0 STEP 1 UNTIL NH DO
02560	BEGIN
02570		NK←N-K;
02580		AA←A[K]+A[NK]; AB←A[K]-A[NK];
02590		BA←B[K]+B[NK]; BB←B[K]-B[NK];
02600		RE←CK*BA+SK*AB;  IM←SK*BA-CK*AB;
02610		B[NK]←IM-BB; B[K]←IM+BB;
02620		A[NK]←AA-RE; A[K]←AA+RE;
02630		DC←R*CK+DC; CK←CK+DC;
02640		DS←R*SK+DS; SK←SK+DS;
02650	END;
02660	END "RTRAN";
02670	
02680	
02690	PROCEDURE RFOUR(REAL ARRAY A,B;INTEGER M,INVERSE);
02700	BEGIN
02710	COMMENT COMPUTES  THE FFT OF 2↑(M+1) REAL DATA POINTS;
02720	INTEGER N,J;
02730	REAL P;
02740	N←2↑M;
02750	IF INVERSE THEN
02760	BEGIN
02770		RTRAN(A,B,N,TRUE);
02780		FOR J←N-1 STEP -1 UNTIL 0 DO
02790		B[J]←-B[J];
02800		FFT(A,B,N,M,N);
02810		FOR J←N-1 STEP -1 UNTIL 0 DO
02820		BEGIN A[J]←0.5*A[J]; B[J]←-0.5*B[J]  END;
02830		REORDER(A,B,N,M,N,TRUE);
02840	END
02850	ELSE
02860	BEGIN
02870		REORDER(A,B,N,M,N,TRUE);
02880		REVFFT(A,B,N,M,1); P←0.5/N;
02890		FOR J←N-1 STEP -1 UNTIL 0 DO
02900		BEGIN A[J]←P*A[J]; B[J]←P*B[J] END;
02910		RTRAN(A,B,N,FALSE);
02920	END;
02930	END "RFOUR";
02940	
02950	
02960	PROCEDURE CFOUR(REAL ARRAY A,B;INTEGER M,INVERSE);
02970	BEGIN
02980	COMMENT  COMPUTES THE FFT OF 2↑M COMPLEX DATA VALUES Xi  IF INVERSE IS TRUE;
02990	INTEGER N,J;
03000	REAL P,Q;
03010	N←2↑M; P←Q←1.0/SQRT(N);
03020	IF INVERSE THEN
03030	BEGIN
03040	Q←-Q;
03050	FOR J←N-1 STEP -1 UNTIL 0 DO B[J]←-B[J];
03060	END;
03070	FFT(A,B,N,M,N); REORDER(A,B,N,M,N,FALSE);
03080	FOR J←N-1 STEP -1 UNTIL 0 DO
03090	BEGIN A[J]←A[J]*P; B[J]←B[J]*Q; END;
03100	END "CFOUR";
03110	
03120	PROCEDURE POWER(REAL ARRAY A,B,C;INTEGER N);
03130	BEGIN
03140	COMMENT THIS COMPUTES THE POWER SPECTRUM OF THE SIN AND COS SERIES IN A,B;
03150	INTEGER I;
03160	FOR I←0 STEP 1 UNTIL N DO
03170	C[I]←SQRT(A[I]↑2 + B[I]↑2);
03180	END "POWER";
03190	
03200	PROCEDURE ARRDIS(REAL ARRAY A; INTEGER N,XPOS,YPOS;STRING ID);
03210	BEGIN
03220	COMMENT DISPLAYS A HISTOGRAM OF THE FIRST N VALUES OF ARRAY A AT XPOS,YPOS;
03230	INTEGER I,J,SP;
03240	INTEGER LY,DY;
03250	REAL MAX;
03260	MAX←0;
03270	FOR I←0 STEP 1 UNTIL N DO
03280	  IF ABS(A[I])>MAX THEN MAX←ABS(A[I]);
03290	MAX←MAX/360;
03300	SP←1024%N;  COMMENT HORIZONTAL SPACING;
03310	AIVECT(XPOS,YPOS); RVECT(1023,0); RIVECT(-1023,0); RVECT(0,360); RIVECT(0,-360);
03320	LY←A[0]/MAX+YPOS;
03330	AIVECT(XPOS,LY);
03340	FOR I←1 STEP 1 UNTIL N DO
03350	BEGIN
03360		DY←A[I]/MAX+YPOS-LY;
03370		LY←LY+DY;
03380		RVECT(SP,DY);
03390	END;
03400	AIVECT(XPOS,YPOS);
03410	FOR I←1 STEP 1 UNTIL 10 DO
03420	BEGIN
03430	  RVECT(0,-15);    COMMENT INSERT HORIZONTAL SCALE;
03440	  RIVECT(26,15);
03450	  RVECT(0,-5);
03460	  RIVECT(26,5);
03470	  RVECT(0,-10);
03480	  RIVECT(26,10);
03490	  RVECT(0,-5);
03500	  RIVECT(26,5);
03510	END;
03520	RVECT(0,-15);
03530	AIVECT(XPOS,YPOS-40);
03540	DPYSST("0      1      2       3      4       5      6       7      8      9      10");
03550	AIVECT(XPOS,YPOS-60);
03560	DPYSST(ID);
03570	END "ARRDIS";
03580	
03590	PROCEDURE DUBDIS(REAL ARRAY A,B; INTEGER N,XPOS,YPOS;STRING ID);
03600	BEGIN
03610	COMMENT DISPLAYS A HISTOGRAM OF THE FIRST N VALUES OF ARRAY A AT XPOS,YPOS;
03620	INTEGER I,J,SP;
03630	INTEGER LY,DY;
03640	MAX←MAX/250;
03650	SP←512%N;  COMMENT HORIZONTAL SPACING;
03660	AIVECT(XPOS,YPOS); RVECT(1023,0); RIVECT(-1023,0); RVECT(0,250); RIVECT(0,-250);
03670	LY←A[0]%10+YPOS;
03680	AIVECT(XPOS,LY);
03690	FOR I←1 STEP 1 UNTIL N DO
03700	BEGIN
03710		DY←A[I]%10+YPOS-LY;
03720		LY←LY+DY;
03730		RVECT(SP,DY);
03740	END;
03750	FOR I←1 STEP 1 UNTIL N DO
03760	BEGIN
03770		DY←B[I]/10+YPOS-LY;
03780		LY←LY+DY;
03790		RVECT(SP,DY);
03800	END;
03810	AIVECT(XPOS,YPOS-40);
03820	DPYSST(ID);
03830	END "DUBDIS";
03840	
03850	REAL PROCEDURE TIMDPY(REAL XPOS,YPOS,RTIM;STRING TFILE);
03860	  BEGIN
03870	  REAL X;
03880	  X←RUNTIM-RTIM;
03890	  IF CHAR1="Y" THEN
03900	    BEGIN
03910	    AIVECT(XPOS,YPOS);
03920	    DPYSST(TFILE&CVS(X));
03930	    END
03940	  ELSE
03950	    OUTSTR(TFILE&CVS(X));
03960	END "TIMDPY";
03970	
03980	INTEGER PROCEDURE FLOW(INTEGER ARRAY DATBUF;INTEGER XPOS,YPOS,M,MAX);
03990	BEGIN
04000	INTEGER ARRAY A[0:512];
04010	INTEGER ARRAY B[0:256];
04020	INTEGER I,J,K,L,BPT,EOF,SEGC,AIFORM,SP,LY,DY;
04030	AIFORM←1;
04040	SP←2;  COMMENT HORIZONTAL SPACING;
04050	  OPEN(AIFORM,"DSK",'10,10,0,0,0,EOF);
04060	  LOOKUP(AIFORM,FILEI&".DMD",0);
04070	  ARRYIN(AIFORM,DATBUF[0],'200);  COMMENT INPUT THE HEADER;
04080	  EOF←0; SEGCNT←0; SEGC←0;
04090	  MAX←2048%MAX;
04100	  WHILE EOF=0 DO
04110	    BEGIN
04120	      ARRYIN(AIFORM,DATBUF[0],DATSIZ);
04130	      IF EOF≠0 THEN
04140		BEGIN
04150		  J←EOF LAND '777777;
04160		  FOR I←J STEP 1 UNTIL 1023 DO DATBUF[I]←0;
04170		END;
04180		BPT←POINT(BPS,DATBUF[0],-1);
04190		FOR I←0 STEP 1 UNTIL 511 DO A[I]←((ILDB(BPT) LSH LEFT)%DI);
04200		FOR I←1 STEP 1 UNTIL 3*DATSIZ%256-2 DO
04210		  BEGIN
04220		    FOR J←0 STEP 1 UNTIL 255 DO  B[J]←((ILDB(BPT) LSH LEFT)%DI);
04230		    FOR K←0 STEP 2↑M UNTIL 255 DO
04240		      BEGIN
04250			DPYSET(DPYBUF);
04260			AIVECT(XPOS+512-K*SP,YPOS); RVECT(511,0); RIVECT(-400,-200);
04270			DPYSST("Segment number "&CVS(J));
04280			LY←A[K]/MAX+YPOS;
04290			AIVECT(XPOS,LY);
04300			FOR L←K+1 STEP 1 UNTIL 511 DO
04310			  BEGIN
04320			    DY←A[L]/MAX+YPOS-LY;
04330			    LY←LY+DY;
04340			    RVECT(SP,DY);
04350			  END;
04360			FOR L←0 STEP 1 UNTIL K DO
04370			  BEGIN
04380			    DY←B[L]/MAX+YPOS-LY;
04390			    LY←LY+DY;
04400			    RVECT(SP,DY);
04410			  END;
04420			DPYOUT(1);
04430		      END;
04440		    FOR J←0 STEP 1 UNTIL 255 DO
04450		      BEGIN
04460			A[J]←A[J+256];
04470			A[J+256]←B[J];
04480		      END;
04490		  END;
04500	    END;
04510	END"FLOW";
     

00010	PROCEDURE LOGDIS(REAL ARRAY A; INTEGER ARRAY D; INTEGER N,XPOS,YPOS;STRING ID);
00020	BEGIN
00030	COMMENT DISPLAYS A LOG FREQUENCY HISTOGRAM OF THE FIRST N VALUES OF ARRAY A AT XPOS,YPOS;
00040	INTEGER I,J,SP;
00050	INTEGER LY,DY;
00060	AIVECT(XPOS,YPOS); RVECT(1023,0); RIVECT(-1023,0); RVECT(0,360); RIVECT(0,-360);
00070	LY←YPOS;
00080	AIVECT(XPOS,LY);
00090	FOR I←2 STEP 1 UNTIL N*210%256 DO
00100	BEGIN
00110		SP←D[I]-D[I-1];	COMMENT HORIZONTAL SPACING;
00120	IF A[I]≤1 THEN DY←YPOS-LY ELSE
00130		DY←100*ALOG10(A[I])+YPOS-LY;
00140		LY←LY+DY;
00150	        IF I=2 THEN RIVECT(SP,DY) ELSE
00160		RVECT(SP,DY);
00170	END;
00180	AIVECT(XPOS,YPOS);
00190	RVECT(0,-15);
00200	RIVECT(92,15);
00210	FOR I←1 STEP 1 UNTIL 7 DO
00220	BEGIN
00230	  RVECT(0,-15);    COMMENT INSERT HORIZONTAL SCALE;
00240	  RIVECT(33,15);
00250	  RVECT(0,-5);
00260	  RIVECT(33,5);
00270	  RVECT(0,-10);
00280	  RIVECT(33,10);
00290	  RVECT(0,-5);
00300	  RIVECT(34,5);
00310	END;
00320	RVECT(0,-15);
00330	AIVECT(XPOS,YPOS-40);
00340	IF SAMINT=2 THEN
00350	DPYSST("0     32       64       128       256      512      1024     2048      4096")
00360	ELSE
00370	DPYSST("0     64       128      256       512      1024     2048     4096      8192");
00380	AIVECT(XPOS,YPOS-60);
00390	DPYSST(ID);
00400	END "LOGDIS";
00410	
00420	PROCEDURE LOGMAG(REAL ARRAY A,B;INTEGER M);
00430	BEGIN
00440	  FOR I←0 STEP 1 UNTIL 2↑M DO
00450	    BEGIN
00460	      A[I]←30*ALOG10(A[I]);
00470	      B[I]←30*ALOG10(A[I]);
00480	    END;
00490	END "LOGMAG";
00500	
00510	 PROCEDURE CEPSTRUM(REAL ARRAY A,B;INTEGER M);
00520	BEGIN
00530	  RFOUR(A,B,M,0);
00540	  LOGMAG(A,B,M);
00550	  RFOUR(A,B,M,1);
00560	END "CEPSTRUM";
00570	
00580	PROCEDURE WINTAB(REAL ARRAY W;INTEGER N);
00590	BEGIN
00600	INTEGER I;
00610	REAL PI;
00620	PI←3.1315926536;
00630	FOR I←0 STEP 1 UNTIL N%2 DO
00640	    W[I]←(1-COS(2*PI*I/N))/2;
00650	END "WINTAB";
00660	
00670	PROCEDURE WINDOW(REAL ARRAY A,B;INTEGER N,K);
00680	BEGIN
00690	COMMENT This procedure applies a Tukey Interim Window to arrays A and B
00700			windowing the k data items from each end ;
00710	INTEGER I;
00720	FOR I←0 STEP 1 UNTIL K DO A[I]←A[I]*W[(N-1)*I/(2*K)];
00730	FOR I←N-1-K STEP 1 UNTIL N-1 DO B[I]←B[I]*W[(N-1)*(N-1-I)/(2*K)]
00740	END "WINDOW";
00750	
     

00010	AIFORM←1; AOFORM←2;
00020	
00030	START:
00040	  IF (TFILEI←STRIN("DATA FILE="))≠NULL THEN FILEI←TFILEI;
00050	  CLOSE(AIFORM);
00060	SAMINT←1;
00070	  IF (TM←ININT("M="))≠0 THEN M←TM;
00080	  IF M>9 THEN GO TO START;
00090	IF M<6 THEN
00100	  BEGIN
00110	    FLOW(DATBUF,-511,0,M,480);
00120	    GO TO START;
00130	  END;
00140	IF M=9 THEN
00150	  BEGIN
00160	    M←M-1;  SAMINT←2;
00170	    OUTSTR("Taken as M=8 and sampling interval of 2"&'15&'12);
00180	  END;
00190	
00200	  N←2↑M; LFILE[1]←DATSIZ%(6*N);
00210	  OPEN(AIFORM,"DSK",'10,10,0,0,0,EOF);
00220	  LOOKUP(AIFORM,FILEI&".DMD",0);
00230	  ARRYIN(AIFORM,LFILE[0],'200);  COMMENT INPUT THE HEADER;
00240	  CHAR1←"Y"; CHAR2←"Y"; CHAR3←"N"; CHAR4←"N"; CHAR5←"N";
00250	  CHAR6←"Y"; LOGFRE←"Y"; CEPS←"Y";
00260	  IF (CHAR7←STRIN("Standard options YorN="))="N" THEN
00270	    BEGIN
00280	      CHAR1←STRIN("Displays   YorN =");
00290	      CHAR2←STRIN("FFT wanted YorN =");
00300	      IF CHAR2="Y" THEN
00310	      LOGFRE←STRIN("Octave scale YorN =");
00320	      CEPS←STRIN("CEPSTRUM YorN =");
00330	      CHAR3←STRIN("Calcom stop YorN=");
00340	      CHAR4←STRIN("Stop anyway YorN=");
00350	      CHAR5←STRIN("Save Spectra YorN =");
00360	      CHAR6←STRIN("Want timings YorN =");
00370	    END;
00380	
00390	IF LOGFRE="Y" THEN
00400	  BEGIN
00410	    FOR I←1 STEP 1 UNTIL 255 DO
00420	        D[I]←1024*ALOG10(I*256/N)/ALOG10(210);
00430	        D[0]←0; D[1]←0;
00440	  END;
00450	
00460	WINTAB(W,N);	COMMENT PREPARE TABLE FOR WINDQW PROCEDURE;
00470	
00480	  IF CHAR5="Y" THEN
00490	    BEGIN
00500	      OPEN(AOFORM,"DSK",'10,0,10,0,0,0);
00510	      ENTER (AOFORM,TFILEI&".POW",0);
00520	    END;
00530	
00540	  EOF←0; SEGC←0; SEGCNT←0;
00550	  SEGTOT←(LFILE[0]*3+2*N-1)%(2*N);
00560	  WHILE EOF=0 DO
00570	    BEGIN
00580	      ARRYIN(AIFORM,DATBUF[0],DATSIZ);
00590	      IF EOF≠0 THEN
00600		BEGIN
00610		  J←EOF LAND '777777;
00620		  FOR I←J STEP 1 UNTIL N-1 DO DATBUF[I]←0;
00630		END;
00640	      BPT←POINT(BPS,DATBUF[0],-1);
00650	      FOR K←1 STEP SAMINT UNTIL 3*DATSIZ%(2*N) DO
00660		BEGIN
00670		  SEGC←SEGC+1;
00680		    DPYSET(DPYBUF);
00690		  IF CHAR1="N" THEN OUTSTR('15&'12&CVS(SEGC))
00700		    ELSE
00710		      BEGIN
00720		        AIVECT(-511,-500);
00730		        DPYSST(DATIME&"   Data file "&FILEI&"   Segment "&CVS(SEGC)&" of "&cvs(segtot)
00740			  &"   SPACING "&CVS(SAMINT)&"   M="&CVS(M));
00750		      END;
00760		  TIMSET;
00770		  UNPACK(A,B,DATBUF,BPT,N,SAMINT);
00780		  TIMDPY(100,-80,RTIM,"  Unpack time = ");
00790		  IF CHAR1="Y" THEN
00800		    BEGIN
00810		      SETFORMAT(2,1);
00820		      DUBDIS(A,B,N-1,-511,256,"Original segment");
00830		      AIVECT(0,420);
00840		      DPYSST("Segment length = "&CVF(N*SAMINT/20));
00850		      IF CHAR2≠"Y" THEN
00860		      DPYOUT(1);
00870		    END;
00880	          IF CHAR3="Y" THEN DPYDMP;
00890		  IF SEGC<SEGCNT THEN
00900		    BEGIN
00910		      WAIT(2);
00920		      IF INCHRS=" " THEN SEGCNT←0;
00930		    END;
00940	          IF CHAR2="Y" THEN IF SEGC≥SEGCNT THEN
00950	            BEGIN
00960		      TIMSET;
00970		      WINDOW(A,B,N,3*N%32);
00980		      TIMDPY(100,-40,RTIM,"WINDOW TIME = ");
00990		      IF CEPS≠"Y" THEN DUBDIS(A,B,N-1,-511,0,"WINDOWED SEGMENT");
01000		      TIMSET;
01010		      RFOUR(A,B,M,0);
01020		      TIMDPY(100,-110,RTIM,"  FFT time = ");
01030		      TIMSET;
01040		      POWER(A,B,C,N);
01050		      TIMDPY(100,-140,RTIM,"  P.Spectrum time = ");
01060		      IF CHAR5="Y" THEN ARRYOUT(AOFORM,C[0],256);
01070		      IF CHAR1="Y" THEN
01080		        BEGIN
01090			  IF LOGFRE="Y" THEN
01100			    LOGDIS(C,D,N,-511,-400,"Power spectrum against log frequency")
01110			  ELSE
01120		            ARRDIS(C,N,-511,-400,"Power spectrum with each major division = "&cvs(1016%samint)&" cycles per second");
01130			IF CEPS="Y" THEN
01140			  BEGIN
01150			    LOGMAG(A,B,M);
01160			    RFOUR(A,B,M,1);
01170			    DUBDIS(A,B,N-1,-511,0,"CEPSTRUM");
01180			  END;
01190			  DPYOUT(1);
01200		          IF CHAR4="Y" THEN
01210		    	    BEGIN
01220		    	      AIVECT(0,-180);
01230		    	      DPYSST("Space bar to continue");
01240			      DPYOUT(1);
01250			     IF(INCHRW)="S" THEN  SEGCNT←CVD(INCHWL);
01260		    	    END;
01270		        END;
01280		      IF CHAR3="Y" THEN DPYDMP;
01290		    END;
01300	        END;
01310	    END;
01320	
01330	  DPYCLR;
01340	  CLOSIN(AIFORM); CLOSO(AOFORM);
01350	  RELEASE(AIFORM); RELEASE(AOFORM);
01360	END "FAST";